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ABSTRACT 

Observations of Lynds Dark Nebula 1221 from the Spitzer Space Telescope 
are presented. These data show three candidate protostars towards L1221, only 
two of which were previously known. The infrared observations also show signa- 
tures of outflowing material, an interpretation which is also supported by radio 
observations with the Very Large Array. In addition, molecular line maps from 
the Five College Radio Astronomy Observatory are shown. 

One-dimensional dust continuum modelling of two of these protostars, IRS1 
and IRS3, is described. These models show two distinctly different protostars 
forming in very similar environments. IRS1 shows a higher luminosity and larger 
inner radius of the envelope than IRS3. The disparity could be caused by a 
difference in age or mass, orientation of outflow cavities, or the impact of a 
binary in the IRS1 core. 

Subject headings: stars: formation, low-mass, L1221 



Introduction 



The story of star formation has l ong included the evolution from a starless core through 
the putative series of collapse stages (jShu et al.lll987l ). These stages include the Class O-III, 
which are defined by various observational signatures and may represent times in t he life 
of the protostar when the envelope and star have certain relative masses ( iLadal 119871 ). The 
Infrared Astronomical Satellite (IRAS) brought about much of what is known about the 
stars forming in dense cores, and, today, the Spitzer Space Telescope reveals new light from 
seen and unseen stars. Through data from Spitzer, new details about the collapse, accretion, 
disk formation, and ultimate birth of stars are being revealed. 

In isolated, low-mass star formation, different cores provide dif ferent laboratories to 
study these processes. Some cores, such as L1014 (jYoung et al.l 12004 ). seem to be forming 
only one star. However, many isolated cores are forming multiple star sys tems. For example 



IRAM 04191+1522 shows two distinct cores with 3 or more protostars (lAndre et al.l 11999 



Dunham et al.ll2006l). and IRAS 16293-2422 is forming a close binary wi th a third star in the 
system (jJorgensen et al.ll2005l ; IChandler et al.ll2005l ; lLoinard et al.ll2007l ). These systems are 
important because they offer views of the impact that other protostars have on the evolution 
of their companions. 

This paper presents new infrared, spectral line, and radio continuum observations of 
Lynds Dark Nebula 1221 (L1221). Descriptions of past observations of L1221 are in §2, 
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including molecular line observations and continuum data from the infrared to millimeter. A 
description of the observations and data reduction of L1221 with the Spitzer Space Telescope 
is in §3. Section 4 has results of the observations including luminosity, classification, and 
other details for the infrared sources detected by Spitzer. Section 5 describes the modelling 
of two protostellar systems in L1221. Finally, §6 offers a discussion and conclusions from the 
modelling. 



2. Lynds Dark Nebula 1221 



Lynd's Dark Nebula 1221 (L1221) was first catalogued by iLvndsl (ll962f ). She listed 
this dark core with an area of 0.020 square degree s and a relative opacit y of 5 (on a scale 
of 1 to 6). The distance to L1221 is 250+50 pc jYonekura et alJll997h . IRAS observed 
an infrared point source towards L1221, IRAS 22266+6845. Sever al authors have rep orted 
on observ a tions of molecular line and (sub)millimeter detections. lYoung et al.l (12006 ) and 
Wu et al.l (120071 ) showed images at wavele ngths of 350-850 /jm , and ICaselli et al.l (120021 ) 
reported their N2H + observations of L 1221. ICaselli et al.l (12002) resolved only o ne extended 
core with a beamsize of 1'. However, lYoung et al.l (120061 ) and IWu et al.l (120071 ) found that 
L1221 has two distinct cores, one in the north (L1221-SMM1) and one towards the south 
(L1221-SMM2), these cores are separated by 50". L1221-SMM1 is at the same position as 
IRAS 22266+6845. 



Both L1221-SMM1 and SMM2 have very similar submillimeter fluxes; lYoung et al. 



( 120061 ) calculated an isothermal mass of 2.6+0.8 M Q for each individual core. These au- 
thors reporte d that L1221-SM M1, which was detected by IRAS, has a larger area than 
L1221-SMM2; IWu et al.l t00l\ ) gave major and minor axes for L1221-SMM1 (46" and 33") 
and L1221-SMM2 (28" and 19"). These axes are measured from the 2-cr contour level. 

Several auth ors have observe d the outflowing; material from LI 221. Based on single- 
dish observations, lUmemoto et al.l (ll99ll ) found a U-shaped outflow, presumably originating 
from IRAS 22266+6845. They concluded that the morphology of the outflow was due to 
interaction with a nearby dense ridge of material but that some external pressure (possi- 
bly magnetic fields) was also required to create the outflow of material seen from IRAS 
22266+6845. iLee et al.l (120021 ) reported, from their interferometric millimeter-wave observa- 
tions (beamsize of 10" x 10"), the J=l— *>0 transition of CO and suggested that the U-shaped 
morphology was due to the interactions of multiple outflows, one of which originates with 
IRAS 22266+6845. Finally. IL~ee~ fc Hoi (j2005h observed L1221 in several molecular lines (CO, 
HCO + , N2H + , and CS) with the Berkeley Illinois Millimeter Array (BIMA) and give 3 mm 
continuum fluxes for MM1 and MM2 (8 and 4 mJy, respectively, with beamsize 4.7" x 3.8"), 
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whose positions are coincident with SMM1 and SMM2. The c ontinuum f l uxes c orrespond 
to masses of 0.02 and 0.01 M for MM1 and MM2, respectively. iLee fc Hoi (bo05h provide a 
detailed model of the kinematics in L1221-MM1. They see evidence of a north-south outflow 
from L1221-MM2, but it is weak and not discussed in detail. They discuss the presence of a 
possible binary system in MM1, which is detected by these Spitzer observations, and drives 
an east-west outflow. The eastern component of this binary is coincident with the position 
of IRAS 22266+6845, but IRAS was unable to resolve the binary pair. ILee fc Hoi (120051 ) 
did not detect any continuum emission around the west object, for which they knew the 
position from th ese Spitzer observations, and concluded there was very little dust around it. 



Lee fc Hoi (120051 ) determined that the observations are best simulated by an infalling, slowly 



rotating, ringlike envelope around the the binary system in L1221-MM1. 



Observations and Data Reduction 



L1221 was observed with the Spitzer Space Telescope (IWerner et al.ll2004f). as part of 



the Legacy Project "From Molecular Cores to Planet-forming Disks" (jEvans et al.l 12003 



c2d). L1221 was observed with two instruments on Spitzer. the Infrared Array Camera 



dFazio et all 120041 . IRAC) at 3.6 (IRAC band 1), 4.5 (IRAC band 2), 5.8 QR AC band 3), 



and 8.0 fim (IRAC band 4) and the Multiband Imaging Photometer for Spitzer ( iRieke et al. 



20041 . MIPS) at 24 (MIPS band 1) and 70 fim (MIPS band 2). The core was observed with 
IRAC on 19 July 2004 (Program ID (PID) 139, AOR key 0005165312) and with MIPS on 
24 September 2004 (PID 139, AOR key 00094287636). 

The IRAC observations of L1221 consisted of 8 pointings arranged in a 2 column by 4 
row grid. Each pointing covers an area of about 5' x 5' for a total area of approximately 
10' x 20'. The area was observed with 4 dithers, or individual images offset by approximately 
10", with exposure times of 12 seconds each giving a total exposure time of 48 seconds. 

The MIPS observation fields were smaller, focusing on the central parts of the LI 221 
core. The 24 fim observations mapped an area of approximately 5' x 15', or a 1 column by 3 
row grid. One cycle (14 frames) of 3 second exposures was used at 24 /im to obtain a total 
exposure time of 42 seconds. The 70 /im camera mapped an area of about 7.5' x 15' in a 
3x3 grid. Each 70 /im pointing had a total exposure time of 90 seconds obtained through 
three 10 frame cycles of 3 second exposures. 

The IRAC and MIPS data were processed through the standard Spitzer Science Center 
pipeline (version S13 for IRAC and MIPS) creating Basic Calibrated Data (BCDs) before 
undergoing further processing in the c2d team internal pipeline. The c2d pipeline improves 
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the BCDs by removing instrumental artifacts in the data. A more complete descr iption of this 



part of the data reduction is available in the c2d data delivery documentation (jEvans et al. 



20071 ). The improved data are then mosaicked to create a single map from the observations 



using the SSC software MOPEX. Sources are extracted and photom etery is performed using 
a modified version of the DoPH OT software (ISchechter et al.lll993l ). The modifications are 
described in lHarvey et al.l (120061 ). 

In addition, the L1221 region was observed with the Very Large Array (VLA) in C 
array configuration on the days of 19 July, 2005 and 25 July, 2005 at 3.6 cm and 6.0 cm. 
Observations were centered on SSTc2d J222807.4+690039, which is called IRS3 in this paper. 
The correlator was setup with four IFs at adjacent frequencies to produce a continuum 
bandwidth of 172 MHz. The data were calibrated and imaged using the standard routines 
of AIPS++. The ob servations and data analysis are similar to the continuum observations 
of L1014 reported in Ishirlev et all ((2003). 

Finally, observations of the N2H + (J = 1 — > 0) and CS (J = 2 — > 1) were taken 
with the 14 meter Five College Radio Astronomy Observatory (FCRAO). The beam size at 
these frequencies is about 55", and the velocity resolution is about 0.07 km/s. The observed 
velocity range for the N2H + is from -5.1 to -3.7 km/s; the line width is 0.68±0.04 km/s. The 
observed velocity range for CS is -5.9 to -2.9 km/s; the line width is 1.5±0.5 km/s. Multiple 
pointings were observed to create maps of the region in the the N2H + and CS transitions, 
which are in th e 3 mm atmospheric w indow. These observations will be described in a 



separate paper (IDe Vries et al 



m prep. 



4. Results 

The images of this region are shown in Figures (H [2], [31 and HI Figure [1] has the IRAC 
3.6, 4.5, and 8.0 /mi observations represented as blue, green, and red. Figure [2] has the 
same color scheme but zoomed in to the region of L1221-SMM1 and SMM2. Three infrared 
sources-IRSl, IRS2, and IRS3-appear to be associated with L1221. IRS1 and IRS2 appear 
unresolved in Figures [1] and [2] because the source is overexposed. The positions of the three 
sources are labeled in Figure [3j which show the IRAC and MIPS data, separately, for the 
3 infrared sources. Finally, Figure HJ has the submillimeter continuum (850 and 350 /mi) 
and FCRAO molecular line observations overlaid on the 8.0 /im greyscale image. Observed 
fluxes for the infrared sources are in Table [H 

SSTc2d J222803.0+690117 (hereafter, L1221-IRS1), was detected at all Spitzer bands 
and was also detected by IRAS (IRAS 22266+6845) and 2MASS (2MASS 22280298+6901166). 
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L1221-IRS1 was also detected in the near- infrared by Hodappl ( 1994 ) and catalogued as an 
outflow driver. This object has [4.5] — [8.0] = 1.4 and [8.0] — [24] = 4.4, whi ch qualifies it 



as a c andidate young stellar object (YSOc) a ccording to th e crite ria given in lHarvey et al. 



( 120061 ); it is also a YSOc after the scheme in lHarvey et al.l (120071 ). Further evid ence of its 
natur e include coincidence with dense gas and evidence that it drives an outflow (ILee fc Ho 
2005|). The fluxes for L1221-IRS1 (1.25 to 850 fim, as in Table ED give L 6o ,=1.8 L and 
T ftn /=250 K, which c lassifies it as a Class I object, according to the scheme developed by 



Myers fc Laddl (119931 ); these wer e calculated using the trapezoidal method of integration. 
The spectral inde x, as defined bvlLadal (119871 ). is a = 0.81, which classifies this object as a 



Class I protostar (IGreene et al.lll994l ). The spectral i ndex is calculated as described in the 
final delivery document for the c2d Legacy Program (jEvans et al.l 120071 ). 



There is a source 7" (1750 AU) west of L1221-IRS1, SSTc2d J22801.8+690119 (hereafter, 
L1221-IRS2). IRS1 and IRS2 were not resolved by IRAS; however, the coordinates for 
the IRAS source are closer to IRS1 than IRS2 (4" versus 17"). The separation between 
IRS1 and IRS2 is about equal t o the resolution of 2MASS, so IRS2 is not included in the 
2MASS catalogs. iHodappI (119941 ) detected several near-infrared sources towards L1221. The 
brightest is clearly associated with L1221-IRS1; they also detect faint emission from IRS2. 
The separation between IRS1 and IRS2 is smaller than the resolution for MIPS at 70 fim, 
17". However, the source position for the 70 /j,m detection is coincident with L1221-IRS1. 
Therefore, an upper limit of about 300 mJy at 70 fim can be assumed for L1221-IRS2; a 
flux larger than this value would cause a shift in the source position towards L1221-IRS2 
(T. Brooke priv. comm.). Additionally, the 24 /xm flux for L1221-IRS2 was not listed in 
the c2d source catalogs, so the IRAC source positions were input to MOPEX, which gave 
the 24 /iin flux of 480 mJy (Brooke priv. comm.). L 1221-IRS2 has [4.5 ] - [8.0] = 2.1 and 
[8.0] - [24] = 2.9, b oth of which qualify it as a YSOc jHarvev et aDl2006h : it is also a YSOc 
after the scheme in lHarvey et al.l (120071 ) . The bolometric luminosity (from 3.6-24 fxm) is 0.4 
L and T M = 450 K (Class I). The spec tral index is a = — 0.05; L1221-IRS2 is considered a 
"flat-spectrum" protostar by this index (IGreene et al.lll994l ). Of course, the spectral energy 
distribution is poorly sampled, and these classifications may change with future observations 
at different wavelengths. 

Spitzer detected a southeastern source at all wavelengths (SSTc2d J22807. 4+690039), 
hereafter L1221-IRS3. This object is about 50" (12500 AU) from L1221-IRS1 and was 
undetected by IRAS or 2MASS. Figures [1] and [2] show 4.5 /im (green) emission emanating 
north and south of IRS3. This emis sion has b e en sh own to indicate outflowing material 
dNoriega-Crespo et al.l 120041) . Indeed, ILee fc Hoi (120051 ) detected a weak, N-S outflow from 



SMM2. 



Lee fc Hoi (I2OO5} ) proposed that this outflow probably emanates from IRS3. 
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L1221-IRS3 has [4.5] - [8.0] = 1. 1 and [8.0] - [ 241 = 5, which qualifies it as a YSOc; 
it is also a YSOc after the scheme in lHarvey et al.l (120071 ). The fluxes give L bo i = 0.8 L 
and Tboi = 68 K over all observed wavelengths (3.6-850 /mi), which classifies it as a Class 
object (<70 K). The spectral index is a = 0.99 (Class I); the spectral index does not define 
the Class 0, so these disparate classifications are not necessarily inconsistent. 

The spectral energy distributions for L1221-IRS1, IRS2, and IRS3 are shown in Figure[5j 
IRS1 shows more emission at shorter wavelengths but is very similar to IRS3 in the longer 
wavelengths, suggesting their envelopes, traced by submillimeter emission, are very similar. 

L1221-IRS3 was the only source of these three that was detected at either 3.6 or 6.0 cm 
with the VLA. This result agrees with th e previous VLA image at 3.6 c m of the L1221 region 
made by IRodriguez fc Reipurthl (119981 ). iRodriguez Reipurthl ( 119981 ) detect an unresolved 
3.6 cm continuum source that is coincident with L1221-IRS3 (22 h 55 m 7.36 s +69° 00' 39.9", 
J2000.0) with a flux of 210 ± 20 /tJy. L1221-IRS3 was detected again in the 2005 VLA 
observations with at 3.6 cm and 6.0 cm with fluxes of 177 ± 24 /tJy and 192 ± 27 /tJy 
respectively. These flux levels are within a factor of 2 of the centimeter fluxes observed 
toward L101 4-IRS, which has 3 .6 cm and 6.0 cm fluxes of 111±8 /tJy and 88±11 /tJy, 
respectivel y (Shirley et al. 12007 ). Also, the 3.6 cm flux for L1221 agrees with the value 
quoted by IRodriguez &: Reipurthl (119981 ). The spectral index between 3.6 cm and 6.0 cm 
is —0.16 ± 0.54. This spectral index is consistent with the index for optically thin free- 
free emission, although a steeper negative or positive spectral index cannot be ruled out. 
Since there is strong evidence from the IRAC images of a molecular outflow from a central 
protostar at L1221-IRS3, the most likely explanation for the observed centimeter emission 
is shock-ionization from i nteraction of the protos t ellar jet from I RS3 with the surr ounding 
material in the envelope (jCuriel et al.lll987L 1198a lAngladalll995l ; Ishang et all 12004 ). 



Finally, there is a fourth infrared source, which was detected from 3.6-70 /mi. SSTc2d 
J222815. 1+685930 is an object about 2.2' southeast from the L1221 core. It is included here 
because it was one of three sources, including L1221-IRS1 and L1221-IRS3, that was detected 
at 70 /mi. The c2d catalogs list t his object as very l ikely a background galaxy based on color - 
magnitude cutoffs from SWIRE fevans et alJl2007l : Irlarvev et al.1120071 : lLonsdale et al.ll2003h . 
Fluxes for this source are given in Table [TJ 



5. 1-Dimensional Models 



Physical models for L1221-IRS1 and L1221-IRS3 were created to best match the fluxes 
from 3.6 to 850 /mi; because the SED of L1221-IRS2 is poorly sampled, there was no at- 
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tempt to model its emission. The models are one- dimensional and included a central star, 
a disk, and envelop e. The disk is included in this one-dimensional model as discussed in 



Dusty (llvezic et al. 



Butner et al.l (11994 ) and Young fc Evans! (120051 ). The radiative transfer was calculated using 



1999). 



The model parameters, which were allowed to vary, include the stellar luminosity (L*), 
the disk luminosity (L D ), the stellar photospheric temperature (T*), the inner radius of the 
envelope (rj), and the power-law density distribution of the envelope (p). 

The other parameters of the model, which were not changed, include the disk inner and 
outer radii (Rj and R ), the power-law surface density distribution of the disk, the power-law 
temperature distribution of the disk (q), properties of dust in the envelope, the outer radius 
of the envelope (r G ), and the total mass of the envelope (M emi ). 

For the disk parameters, default values for all except the disk luminosity, Ld, have been 
adopted. The surface density power law index is p = 1.5, and the temperature power law 
index is q = 0.5. The for mer is in accordance with the density structure for a disk in vertical 
hydrostatic equilibrium (IChiang fc Goldreichlll9970 ; and the latter is chosen to simul ate the 
effects of flaring and disk accretion (IButner et al.l Il994l ; iKenyon &: Hartmannl 119871 ). The 
inner radius of the disk is set to where the dust is heated to greater than 2000 K and is, 
presumably, destroyed. The outer radius is assumed to be 5 AU; this is similar to the dis k 
radius, based on the centrifugal radius, used for IRAM 04191+1522 (IDunham et al.ll2006l ). 
The effect of the disk outer radius (R Q ) is discussed in §5.31 The scale of the power-law 
density distribution was set by assuming th e disk has a mass of 0.005 M , which is similar 
to the assumed disk for IRAM 04191+1522 (IDunham et al.ll2006l ). a source that seems to be 
similar to L1221-IRS3. In fact, the mass of the disk has little effect on the model. The disk 
luminosity, Lp, arises because of accretion onto the disk and is a free parameter in these 
models. Bas ically, Lb has to be tu ned in order to match the 24 and 70 /im observations. As 
discussed in lYoung fc Evanj (120051 ). Lb is not indicative of the mass of the disk; it is simply 
set to a particular value so as to match the observations. 

The 850 fxm flux is least susceptible to geometric effects because the envelope is optically 
thin at this wavelength. Therefore, this flux was used to constrain the mass of the envelope. 
The fiducial density was set, depending on the value for the inner radius, such that the 
modelled 850 ^m flux would match the observed flux. 

Because the shape of the envelope's density distribution is not known, a power-law 
density distribution, n(r) = (j r;^ ( r f = 1000 AU), with indices of p = 1.5 or p = 2.0 



is assumed. lYoung et al.l (120031 ) found a median p of 1. 8; indeed, most studies have found 
an average power-law index to be between 1.5 and 2.0 (iShirley et al.ll2002l ; iMotte fc Andre 
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200l|). 



The outer radius of the envelope is chosen to be r Q = 5000 AU. The outer radius is 



not easily constrained from submillimeter observations. However, IWu et al.l (120071 ) gave the 
major axes of IRS1 and IRS3 (46" and 28"; 12000 and 7000 AU) at the 2-a level. These 
sizes correspond to outer radii of 6000 and 3500 AU. There f ore, t he outer radius is set to the 
average of these values (5000 AU) for both cores. IWu et al.l (120071 ) also warned that SHARC- 
II is not sensitive to extended emission, so the reported values potentially underestimate the 
submillimeter size of the cores. However, the cores are separated by 50" in the plane of the 
sky, which corresponds to 12500 AU. Then, it is reasonable to assume the outer radii do not 
exceed 6000 AU. 



For the dust properties, the opacities calculated by lOssenkopf Hennind (119941 ) from 
the fifth column of Table 1 in their paper (so-called OH5 dust) were used. The s e dat a have 
been extended to a greater range of wavelengths as described in 
scattering of light by the dust grains is ignored, as discussed in 
impact of this assumption is discussed in §5.31 



Young & Evaraj(|2005() . The 



Young & Evansl (1200a ). The 



in 



The envelope is heat ed externally by the interstellar radiation field (ISRF) . As described 



Young fc Evansl (120051 ). the ISRF has been attenuated by the opacity due to lDraine &: Lee 
(119841 ) dust with Ay = 3. This simulates the effect of low density material in the surrounding 



environs of a star-forming core. 

To determine the best-fit model, the reduced x 2 was calculated over all Spitzer wave- 
lengths as such 



x 2 



1 ^ [s° bs (\ 

k 2-*> 



3.6 < Aj < 70fim 



MA.)] 2 

The degrees of freedom (k = n — m) is the difference between the number of free parameters 
(m = 4) and the number of data points in the SED (n = 9). 

Then, four different grids were examined to find the parameters that give the lowest x 2 
values. These four grids cover the following parameter spaces: — Ld> — r^, — 7"j, and 
T* — L*, for grids 1 through 4, respectively. 

The x. 2 values for grids 1 and 2 are plotted as greyscale in Figure [6] for L1221-IRS1 
and Figure [7] for L1221-IRS3. Grids 3 and 4, which are not shown, pinpoint values for 
Ti and that are similar to those found in grids 1 and 2. In addition, grids 3 and 4 
probe the effective temperature of the protostar. These models are insensitive to stellar 
temperature. The envelope is optically thick at the shorter wavelengths, so all of the stellar 
photospheric radiation is reprocessed. The stellar temperature, then, is set at 3000 K, which 
is representative for a late-M dwarf. Because the x 2 plots for Grids 3 and 4 are redundant, 
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they are not shown, but the minimum x 2 values for several parameters from all grids are 
plotted in Figures [8] and [91 



5.1. L1221-IRS1 

L1221-IRS1 was detected by IRAS, Spitzer, 2MASS, SHARC-II, and SCUBA. How- 
ever, the near-infrared data are not well fitted by a one-dimensional model, so the 2MASS 
observations are not used in calculating x 2 (as given by equation 1). 

L1221-IRS1 has a nearby companion, L1221-IRS2, about 7" away, which corresponds to 
1750 AU. However, L1221-IRS2 is much less luminous than L1221-IRS1; it was not resolved 



at 70 ^m, but there is an upper limit on the 70 jum flux of 300 mJy Additionally, iLee fc Ho 



( 120051 ) detected no continuum 3 mm emission from L1221-IRS2 an d concluded that there 



was probably very little dust surrounding this object. IHodappI ( 119941 ) detected emission from 



both IRS1 and IRS2, though the emission from IRS2 was very faint. 

Because IRS2 is much less luminous, it is not included in the model of IRS1 and its 
surrounding envelope. It is, of course, possible that L1221-IRS2 is not even associated with 
the core and is simply a background object. However, it is a YSOc, so its impact as a 
potential binary companion will be discussed in section 6. 

Initially, values for the model parameters that produced a reasonably good fit to the 
data were assumed: = 1000 AU, L* = 1 L Q , Ld = 1 L , and T* = 3000 K.; these are 
the default values for the free parameters. When allowed to vary, these parameters were set 
in the following ranges: r, = 250 — > 2500 AU, L* = 0.5 — > 5 L Q , L D = 0.5 — > 5 L , and 
T* = 1000 20000 K. 

The fiducial density, at r/ = 1000 AU, was set such that the total envelope mass was 
always a certain value to match the 850 fim flux (2 Jy). For the models with p = 1.5 , the 
mass is 1.1 M Q ; when p = 2.0 , M env = 0.9 M Q . Also, when the envelope has different inner 
radii, the fiducial densities must be adjusted to produce these envelope masses. For inner 
radii from 250 to 2500 AU, the range of fiducial densities were 1.8 x 10 6 to 2.7 x 10 6 cm -3 
for the p = 1.5 models and 2.2 x 10 6 to 4.2 x 10 6 cm -3 for the p = 2.0 models. 

Figure E] shows greyscale plots with the x 2 values of models with p = 1.5 and p = 2.0. 
Grid 1 shows a range of appropriate luminosities for the star and disk. Grid 2 effectively 
pinpoints one valid inner radius and stellar luminosity. 

For each set of models, the minimum x 2 values for the total internal luminosity (L int ) 
and rj are plotted in Figure [8J These graphs show which of these values gave the best x 2 
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and were used in determining the best-fit parameters. 

To determine the appropriate model parameters, the inner radius is selected from Grid 
2. Then, the inner radius, which gives the lowest x 2 , is taken from Grid 3 (see Figure [S]). 
Averaging the values for ri gives the best-fit model parameter as reported in Table [2j There 
is some ambiguity in the knowledge of L* and Lp,. As shown in Grid 1 of Figure El the best- 
fit stellar luminosity depends upon the chosen disk luminosity. Therefore, the sum L* + Lp> 
(L int ) is constrained with Grid 1 and equally apportioned for and Lp>. 

The best-fit parameters are given in Table [2j The model with p = 2.0 has a slightly 
better y 2 value; however, both models have fairly high x 2 values (~50). For both p — 1.5 
and p = 2.0, the inner radius (r^) is 1000 AU; the total internal luminosity is 2.6 L . The 
spectral energy distributions for these two models are shown in the upper panels of Figure [TU1 
The dotted line in these plots represents the spectral energy distribution (SED) of the star 
and disk, excluding the envelope; the solid line is the SED of the star, disk, and envelope 
system. The error bars represent the observed fluxes. 

The models often did not match the MIPS data, so X^mipsi which only includes Aj =24 
and 70 /xm, is used to find the model that best fits the MIPS data. L int was held constant 
while Lp> was varied to optimize Xmips w hile the envelope parameters (r i; p, and rif) are 
the same as in Table [2j These models are represented by the long dashed line in Figure [TOl 
For p = 1.5, = 1.3 and Lp> = 1.9 L ; for p = 2.0, = 1.3 and = 1.7 L ; the best 
Xmips varues were for the p — 1.5 model. In both cases, the 70 /im observations are better 
matched, but the 24 /im is still underestimated by the model. As Lp> increases, the 70 jum 
flux increases at a greater rate than the 24 /im flux because of the higher extinction at 24 
/iiii. Therefore, as Lp> is increased, the 70 /mi flux is increased, but the 24 /im remains about 
the same. Of course, the IRAC fluxes are also overestimated by these models, since their 
values were not considered in calculating Xmips- This effect is prevalent in all four models. 
A more realistic and two-dimensional treatment of the disk is necessary to remove this effect. 



5.2. L1221-IRS3 

L1221-IRS3 is the only infrared source detected in the southern submillimeter core, 
L1221-SMM2. This is the first report of its infrared detection. Modelling of IRS3 follows 
closely the process for IRS1. The same default values for disk parameters (except L D and 
Ri), dust opacity (OH5), power-law density profiles (p = 1.5 and 2.0), and heating by the 
ISRF as used for IRS1 were kept for IRS3. 

Initially, values for the model parameters that produced a reasonably good fit to the 
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data were assumed: r{ = 100 AU, L* = 0.4 L Q , Lp = 0.4 L , and T* = 3000 K; these are 
the default values for the free parameters. When allowed to vary, these parameters were 
set in the following ranges: = 50 — > 750 AU, = 0.1 — > 2 L Q , = 0.1 — > 2 L Q , and 
L* = 500 -> 5000 K. 

Then, the fiducial density (n/) was set to match the modelled and observed 850 jum 
fluxes. The 850 /im flux is dependent on the dust mass and temperature. The models with 
different values for p have differing temperature distributions. Therefore, one must choose 
2 envelope masses that create an 850 //m flux that best matches the observations. For the 
p = 1.5 models, the envelope mass is 1.75 M . When p = 2.0, M env = 1.25 M Q . The steeper 
density profile has higher temperatures, so a smaller dust mass is required to produce the 
850 /im flux. 

Because the inner radius changes in these models, different fiducial densities are required 
to create an envelope with these masses. For inner radii from 50 to 750 AU, the range of 
fiducial densities (at 77 = 1000 AU) were 2.8 x 10 6 to 2.9 x 10 6 cm 3 for the p = 1.5 models 
and 3.0 x 10 6 to 3.5 x 10 6 cm 3 for the p = 2.0 models. 

The x 2 two-dimensional plots are in Figure [71 The plots of minimum x 2 are in Figure [9j 
Analyzing these data in the same way as for IRS1, the best- fit parameters were found, which 
are given in Table [21 are as follows. For p = 1.5 , = 100 AU, and Li nt = 0.4 L Q . For 
p = 2.0 , rj = 150 AU and L int = 2.2 L . The model with p = 2.0 has a slightly better x 2 
value, though both models have fairly high x 2 values. The spectral energy distributions for 
these two models are shown in the lower panels of Figure [101 

As discussed for IRS1, models considering just Xmips were calculated and are shown 
by the long-dashed lines in Figure [TUl The best-fit models have the same parameters for the 
envelope (r*j, p, and n/) as those in Table [21 but their disk luminosities are varied. The new 
parameters are as follows: for p = 1.5, L int =0.2 and Lo =1.0 L ; for p = 2.0, L int =1.1 and 
Ld =0.4 L . The best Xmips va hies were for the p = 1.5 model. 

The modelled L int is significantly different than the observed bolometric luminosity (0.8 
L ). This discrepancy highlights the importance of far-infrared observations. For the p = 2.0 
model, the luminosity between 75 and 350 /im is 1.9 L . However, this spectral range, which 
is the peak of the SED, is not sampled with observations, so the modeled L« is much larger 
than the observed L bol . 
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5.3. Constraint of Internal Luminosity 



These models are one-dimensional and insufficient to fully describe the physical nature 
of these systems. One-dimensional models are not equipped to simulate the effects of a 
more realistic disk or scattering by dust grains. These shortcomings cause an uncertainty in 
the modelling of fluxes at shorter wavelengths, especially in the near- infrared. As such, the 
NIR fluxes have not been included in the x 2 calculations. However, one-dimensional models 
are useful because they constrain the luminosity of the internal source. Either including 
scattering of light or altering disk parameters does not affect the conclusions about internal 
luminosity. 

Because dust grains preferentially forw ard scatter light, the on e -dimensional code, Dusty , 
is unable to properly account for scattering (jYoung &: Evansll2005l ) . iDunham et al.l (jin prep.1 ) 
have created dust opacities that account for scattering by adding t he scattering and absorp - 
tion opacities. This is not an entirely reasonable assumption, but IDunham et al.l (jin prep.1 ) 
have shown that this method offers good agreement with two-dimensional models, which are 
able to treat scattering in a more correct manner. 

Figure [11] shows the best-fit models, whose parameters are given in Table [2j as a dotted 
line. The s olid line in Figure ITT1 sh ows the best-fit model when these new opacities are 
used (as in IDunham et al.l (jin prep.1 )). The new opacities allow for a better fit at shorter 



wavelengths, especially in the near- infrared. However, the internal luminosity for IRS1 and 
IRS3 are changed very little; the inner radii for both sources are increased by, at least, a 
factor of 2. 



Dunham et al.l (120061 ) concluded, similarly, that the luminosity of the internal source 
is well-constrained by one- dimensional models. The results for L int of IRS1 and IRS3 are 
robust, but the inner radii are dependent on a variety of parameters related to the disk, dust 
opacities, and other two-dimensional effects. 

In addition, the disk is important to the overall fit of the model, but certain parameters 
of the disk, such as its outer radius, do not have a large impact on the best-fit models. A 
small disk outer radius (5 AU) has been adopted for the models of IRS3 and IRS1, but the 
radius of the disk has little effect on the modelled SED. Figure [12] has the best-fit models 
for IRS1 and IRS3 (with p — 1.5). In addition, this figure shows models with disk outer 
radii of 50 and 500 AU. Though the disk spectra (dashed line) are markedy different, the 
SED of the protostar, disk, and envelope (solid line) are all about the same. Clearly, the 
disk plays an important role in the mid-infrared portion of the SED; in fact, these models 
underestimate the 24 /im flux. However, further observations to constrain disk parameters 
and higher-dimensional models are needed for more certain conclusions. Observations with 
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Sofia and Herschel will provide better sampling of the SED, and submillimeter interferometry 
will be useful to constrain the disk masses and sizes in this system. 



In comparison, iDunham et al.l (120081 ) found a relationship between L int and the flux 
at 70 (Ltm; the relationship is found to be reliable within a factor of ~ 2. Based on the 70 
Atm flux for IR S 1, the derived L int is 1.1 L Q , while the model suggests L int = 2.6 L Q . The 
Dunham et al.l (120081 1 relationship derives 0.8 L for IRS3, while the model gives a range 
for L int from 0.4 to 2.2 Lm, dep e nding on the density profile for the envelope. Given the 
uncertainty in the lDunham et al.l (120081 ) relationship, the modeled luminosities are consistent 
with those derived from the 70 /im flux. 

In conclusion, results for the internal luminosity are fairly robust. They are not affected 
by changes either in the disk or the effects due to scattering, even though these changes can 
drastically affect the near-infrared portion of the overall spectrum. 



Discussion and Conclusions 



The L1221 star-forming region presents an example of stars forming differently in very 
similar environments. L1221-SMM1 and SMM2 are, presumably, both part of the same 
region and are separated, in the plane of the sky, by approximately 50". Their 850 /im fluxes 
are identical, so their total masses are also similar, though their submillimeter sizes and, 
possibly, density profiles, are different. Also, the two cores ap pear to be a part of the same 
molecular core, as in Figure HJ and dark cloud (jLynda Il962l ) . These two cores, because of 
their close proximity (50" = 12500 AU), are likely to be affected by the same large-scale 
dynamical effects due to nearby stars and astronomical events; they are also likely to have 
been born out of a very similar makeup of materials. In short, L1221-SMM1 and SMM2 
should provide similar environments for the formation of stars. 

L1221 is not unique. A number of c ores have been f ound to be home to multiple pro- 
tostars in different stages of formation. iDuchene et al.l (120071 ) found multiplicity rates of 
30-50% among protostellar systems with separations of up to 1400 AU; the protostars in- 
cluded Class I and flat-spectrum sources in several mo l ecular clouds. IRAM 04191+1422 has 
several similarities with L1221, as well . lAndre et al.l (119991 ) discovered this very-low lumi- 
nosity object (VeLLO) JPunham et al1l2006h about r (8400 AU) from IRAS 04191+1523, 
a Class I protostar. Like L1221-IRS1, IRAS 04191+1523 also has two sources, unresolved 
by IRAS, that ar e 6.5" apart. These sources were resolved by IRAC and 2MASS but not 
MIPS or SCUBA (IDunham et al.l 120061 ). This star-forming region is similar to L1221 in that 
it appears to have several protostars at different stages of evolution. L1251B is another 
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star- forming region th at is hom e to different cla s ses of protostars including four Class O/I 
stars (ILee et al.ll2006l ). Finally, lYun &: Clemens! (119951 ) found two near- infrared sources in 
the star-forming core CB230. T hese supposed proto stars are sep arated by 11". O ther ex- 
amples of pairs include BHR 71 faourke et al.ll200lh and CG30 Jchen et al.lhoosh . These 
cores are examples of sources similar to the L1221 star-forming region in that they are all 
forming multiple protostars that are, presumably, at differing stages of evolution. 

IRS1 and IRS3, the dominant infrared sources in SMM1 and SMM2, are clearly not 
as similar as their host cores. First, their luminosities are different. L1221-IRS1 has an 
internal luminosity, including the disk and stellar components, of 2.6 L , while L1221-IRS3 
has Li n t = 2.2 or 0.4 L , depending on the density profile. The observed infrared luminosity 
(with A < 70 fim) for IRS 1 is 1.6 L and for IRS3 is 0.4 L . 

The luminosity is the result of accretion onto the disk and star. Since each core likely 
has similar temperature and turbulence, the accretion rates for IRS1 and IRS3 are probably 
very similar. Because accretion luminosity is directly proportional to the protostar's mass, 
IRS1 is potentially about 4 times more massive than IRS3. 

However, IRS3, for the two best-fit models, has a range of possible luminosities from 
0.4 to 2.2 L . Additional observations in the wavelength range from 70 to 350 fim (such as 
with Herschel) would easily discern the appropriate internal luminosity since the models are 
quite different in this portion of the SED. 

Current submillimeter (450 and 850 fim) maps, made with the scan map technique on 
SCUBA, have poor signal-to-noise ratios and are insufficient to determine the shape of the 
density profile for these cores. More sensitive observations of the e xtended submillime ter 
emission are needed to determine the density profile of the envelope (jShirley et al.ll2002l ). 



The envelopes' inner radii of the L1221-IRS3 and L1221-IRS1 models are distinctly 
different. L1221-IRS3 requires = 100 or 150 AU; L1221-IRS1 needs a much larger inner 
radius (1000 AU). Two scenarios can explain this disparity. First, L1221-IRS3 might be at 
a much earlier stage in its colla pse, and the inner radius has not expanded as predicted by 



Terebey. Shu. &: Cassen I (jl_984j). Another alternative is that the binary possibly in L1221- 
SMM1 (IRS1 and IRS2) could have evacuated the cavity of some f raction of the m a terial . 
The separation of this binary is 7", which corresponds to 1750 AU. IJeirgensen et al.l (120051 ) 
found a similarly large inner radius for IRAS 16293-2422 and repor t ed th at the radius of the 
inner cavity was comparable to the centrifugal radius. ILee &: Hoi (120051 ) reported a "hole" 
in the N 2 H + emission, which could suggest a cleared out c avity in the cent e r of t he core. 
However, th is could also be a chemical effect as discussed in IJeirgensen et al.l (120041 ) and as 



observed by I J0rgensenl (120041 ) . 
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For IRAS 16293-2422, Ijprgensen et all (120051 ) suggested a two-dimensional model that 



used a much smaller inner radius when outflow cavities were instituted and when the pro- 
tostar was viewed down the outflow cavity. They found that the mid-infrared fluxes were 
greatly increased when observed down the outflow cavity A similar scenario might apply 
for L1221 IRS1 and IRS3. IRS1 might be observed pol e-on while IRS3 is edge-on. However, 



because of the bipolar nature of IRS1 (ILee fc Holl2005l ). this is unlikely. 



The chief results of these models include these conclusions about the luminosities and 
envelope inner radii. The luminosity of each core is easily constrained with a one- dimensional 
model because most of the lumonisity is emitted in the far-infrared, which is largely unaf- 
fected by the central star's near-infrared spectrum; the one- dimensional model is sufficient 
also because the re-radiated emission is thin. Given the same data, even a more complex 
two-dimensional model should reach the same conclusions regarding the luminosity of IRS1 
and IRS3. Prediction of the envelope's inner radii, on the other hand, is a less robust conclu- 
sion. The central protostar's spectrum does have a large impact on the inner radius required 
to match the data. However, the models do suggest that the inner radius of IRS1 is much 
larger than the inner radius of IRS3. These results still hold even when the disk is altered 
or scattering is included in the model. 

Finally, IRS1 and IRS2, the potentially binary companion to IRS1, show different SEDs. 
IRS1 is detected at submillimeter wavelengths, while IRS2 shows no clear evidence of long- 
wavelength emission. Possibly, IRS2 is more evolved than IRS1 and has stopped accreting 
material from the surrounding envelope. Perhaps, IRS1 and IRS2 exhibit such different 
properties because IRS1 is dominating the accretion. In fact, this might be common in 
binaries that are close enough; initially, one gets bigger and then dominates the accretion 
for the duration of their formation. 
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Fig. 1. — Three-color IRAC image of L1221. Red corresponds to 8.0 fxm, green is 4.5 /im, and 
blue is 3.6 /im. The sources are labeled in Figured IRS1/IRS2 is the brightest source; they 
appear unresolved because the source is overexposed in this image. IRS3 is the small, green 
object southeast of IRS1/2. An arc of emission extends from L1221-IRS1 to the northeast. 
L1221-IRS3 shows a signature of a north-south outflow. 
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Fig. 2. — Three-color IRAC image of L1221. The color scheme is the same as in Figured] 
the image is zoomed in to show detail of the inner region. IRS1 and IRS2 are in the brightest 
source; they appear unresolved because the source is overexposed in this image. IRS3 is the 
green object southeast of IRS 1/2. 
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Fig. 3. — The IRAC and MIPS images are shown for L1221. The maps are shown with 
logarithmic scaling to show IRS1, IRS2, and IRS3. 
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Fig. 4. — Millimeter observations are overlaid, as contours, on the 8 /im greyscale image. 
The upper left panel shows 350 yum, upper right is 850 /mi, the lower left panel has CS 
(J = 2 — ► 1), and the lower right panel has N 2 H + ( J = 1 — > 0). 
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Fig. 5.— Spectral energy distributions for L1221-IRS1, IRS2, and IRS3. IRS1 and IRS3, 
shown as solid lines, were modelled while IRS2 (dashed line) was not. 
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Fig. 6. — The x 2 values for models of IRS1 with p = 1.5 and p = 2.0. White areas denote 
lower x 2 values. The default values, except when allowed to vary, are L* = 1.0 L , L D = 1.0 
L , n = 1000 AU, and T* = 3000 K. 
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Fig. 7. — The x 2 values for models of IRS3 with p — 1.5 and p = 2.0. White areas denote 

lower x 2 values. The default values, except when allowed to vary, are L* = 0.4 L , L D = 0.4 

L , n = 100 AU, and T* = 3000 K. 
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Fig. 8. — The minimum x 2 values plotted from the data shown in Figure [6j 
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Fig. 9. — The minimum x 2 values plotted from the data shown in Figures 
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Fig. 10. — The best-fit models for differing power-law indices have different parameters as 
shown in this figure for both L1221-IRS1 and L1221-IRS3. Best-fit parameters are in TableEJ 
The dotted line represents the SED of the star and disk; the solid line is the SED of the star, 
disk, and envelope. The dashed line is the model that best-fits only the MIPS observations. 
The error bars are the observed fluxes. 
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Fig. 11. — The effects of scattering have been included in the SEDs shown in this figure. 
The dotted line shows the models from Figure [10] while the solid line shows best-fit models 
when scattering is included. The parameters for these best-fit models are as labeled here. 
The internal luminosity is mostly independent of the inclusion of scattering; the inner radii 
increase by a factor of 2 or more when scattering is included. 
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Fig. 12. — SEDs for IRS1 and IRS3 (p = 1.5) in which the disk outer radius has been 
changed. The blue shows the SED with R = 5 AU, red is R Q = 50 AU, and green is 
R = 500 AU. Black error bars are the observed data. 



-32 - 



Table 1. Spectral Energy Distribution 



Source Name 


Lambda 


F„ 


Aperture 




(im 


mjy 


arcsec 



SSlcid J222803. 0+690117 


1.25 


0.374+0.037 


2.5 


(L1221-lKbl) 


1.65 


3.91+0.39 


2.5 


T Obs I Q T 


n 17 


I T V 

I I .ortl. / 


z.o 


1 bo! - 250 K 


o.O 


1 nn_l_i i 


1 To 
1. 1 


a = 0.81 


/I c 

4.0 


1 /t n_l_i c: 
14y±lo 


1 TO 
1. 1 




O.O 


Qno_l_Qn 


1 no 
1.9 




8.0 


38/ ±39 


2.0 a 




24 


i o /i n_l_ i o 7 
1940±187 


5.7 a 




70 


6940+641 


17 a 




350 


9300+1400 


40 




450 


10400+6500 


40 




850 


2000+200 


40 


SSTc2d J222801. 8+690119 


3.6 


83.6+ 8.4 


1.7° 


(L1221-IRS2) 


4.5 


137+ 14 


1.7 a 


L^f=0.4 L 


5.8 


324+32 


1.9 a 


Tg o f=450 K 


8.0 


381+38 


2.0 a 


a = -0.05 


24 


368+35 


5.7 a 




70 


< 300 


17 a 


SSTc2d J222807.4+690039 


3.6 


0.567+0.060 


1.7 a 


(L1221-IRS3) 


4.5 


3.1 ±0.3 


1.7 a 




5.8 


4.96+0.50 


1.9 a 


Tgjf =68 K 


8.0 


3.84+0.4 


2.0 a 


a = 0.99 


24 


47.5+4.4 


5.7 a 




70 


5080+469 


17 a 




350 


6400+1000 


40 




450 


11300+6100 


40 




850 


2000+300 


40 


SSTc2d J222815. 1+685930 


3.6 


0.208+0.02 


1.7° 




4.5 


0.284+0.03 


1.7 a 


Kot =0-003 L 


5.8 


0.38+0.04 


1.9 a 


T£?=187 K 


8.0 


0.613+0.06 


2.0 a 


a = 0.25 


24 


2.13+0.27 


5.7 a 




70 


24.1+3.8 


17 a 



FWHM of the Spitzer point-spread profile. 
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Table 2. Model Parameters 





Parameter 


Accepted Value 


L1221-IRS1 


V 


1.5 






5xl0 6 cm" 3 






1000 AU 




Lint 


2.6 Lq 






3000 K 


L1221-IRS1 


p 


2.0 




n f 


7xl0 6 cm- 3 




n 


1000 AU 




Lint 


2.6 Lq 






3000 K 


L1221-IRS3 


p 


1.5 




n f 


5xl0 6 cm- 3 




n 


100 AU 




Lint 


0.4 Lq 






3000 K 


L1221-IRS3 


p 


2.0 




n f 


7xl0 6 cm^ 3 




Ti 


150 AU 




Lint 


2.2 Lq 




T« 


3000 K 



